{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diffraction Limited PSF and MTF\n",
    "\n",
    "In this example, we will show how to calculate the diffraction limited PSF and MTF for a circular aperture.  We will also compare the numerically derived results from `prysm` with the analytical forms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from prysm import Pupil, PSF, MTF\n",
    "from prysm.psf import airydisk\n",
    "from prysm.otf import diffraction_limited_mtf\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.style.use('bmh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# system parameters\n",
    "epd = 1\n",
    "efl = 10\n",
    "fno = efl / epd\n",
    "wavelength = 0.5\n",
    "\n",
    "p = Pupil(wavelength=wavelength, dia=epd, samples=256)\n",
    "psf = PSF.from_pupil(p, efl)\n",
    "u, sx = psf.slices().x\n",
    "\n",
    "# calculate the analytical version, and the difference between the two\n",
    "analytical_psf = airydisk(u, fno, wavelength)\n",
    "psferr = (analytical_psf - sx)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10,4))\n",
    "psf.slices().plot('x', fig=fig, ax=ax1, xlim=50)\n",
    "ax1.plot(u, analytical_psf, ls=':', c='k', label='Analytic', zorder=3)\n",
    "ax1.legend()\n",
    "ax2.plot(u, psferr, lw=3)\n",
    "ax2.set(xlim=(-50,50), xlabel=r'Image Plane X [$\\mu m$]', ylabel='Intensity Difference an - numerical')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can see that the differences manifest below the fourth decimal place.  It might be interesting to see the RMS error,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from prysm.util import rms\n",
    "\n",
    "rms(psferr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This error is over a 1D slice.  If calculated over the whole image plane it would be much smaller.\n",
    "\n",
    "Next, we consider the MTF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the MTF and its error\n",
    "mtf = MTF.from_psf(psf)\n",
    "nu, m = mtf.slices().x\n",
    "m_analytic = diffraction_limited_mtf(fno, wavelength, frequencies=nu)\n",
    "mtferr = (m_analytic - m)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10,4))\n",
    "mtf.slices().plot('x', fig=fig, ax=ax1)\n",
    "ax1.plot(nu, m_analytic, ls=':', c='k', label='Analytic', zorder=3)\n",
    "ax1.legend()\n",
    "ax2.plot(nu, mtferr, lw=3)\n",
    "ax2.set(xlim=(0,200), ylim=(-0.0005,0.0005), xlabel='Spatial Frequency [cy/mm]', ylabel='MTF Difference [a.u.]')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again the error is at the fourth decimal place.  The RMS may be interesting again,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rms(mtferr)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}